Existence of short-time approximations of any polynomial order for the computation 

of density matrices by path integral methods* 
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In this article, I provide significant mathematical evidence in support of the existence of di- 
rect short-time approximations of any polynomial order for the computation of density matrices 
of physical systems described by arbitrarily smooth and bounded from below potentials. While 
for Theorem [5] which is "experimental," I only provide a "physicist's" proof, I believe the present 
development is mathematically sound. As a verification, I explicitly construct two short-time ap- 
proximations to the density matrix having convergence orders 3 and 4, respectively. Furthermore, 
in Appendix B, I derive the convergence constant for the trapezoidal Trotter path integral tech- 
nique. The convergence orders and constants are then verified by numerical simulations. While the 
two short-time approximations constructed are of sure interest to physicists and chemists involved 
in Monte Carlo path integral simulations, the present article is also aimed at the mathematical 
community, who might find the results interesting and worth exploring. I conclude the paper by 
discussing the implications of the present findings with respect to the solvability of the dynamical 
sign problem appearing in real-time Feynman path integral simulations. 
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I. INTRODUCTION 

In the path integral formulation, the density matrix 
of a thermodynamic system is expressed as the expected 
value of a functional of the Brownian motion by means 
of the Feynman-Kac formulaii2iii 



p(x,x';(3) 



Eexpl -(3 / V 



x r (u) + o-B[ 



du 



Pf p (x,x';[3) 

(1) 

Here, p(x,x';(3) is the density matrix for a onc- 
dimensional canonical system characterized by the in- 
verse temperature (3 = l/(fcsT) and made up of identical 
particles of mass mo moving in the potential V{x). The 
stochastic element that appears in Eq. iTI}, {B®, < u < 
1}, is a so-called standard Brownian bridge, defined as 
follows: if {B u , u > 0} is a standard Brownian motion 
starting at zero, then the Brownian bridge is the stochas- 
tic process {B u , < u < 1| B\ = 0}, i.e., a Brownian mo- 
tion conditioned on the event B\ = 0. A Brownian bridge 
can be realized as the process {B u — uBi,Q < u < l}. 4 
For additional information on Brownian motion and its 
relation to the Feynman-Kac formula, the reader is ad- 
vised to consult Appendix A as well as the cited bib- 
liography. To complete the description of Eq (Tj^), we 
set x r (u) — x + (x 1 — x)u (called the reference path), 
o~ = (h 2 [3/mo) 1 / 2 , and let Pf p (x, x'; /?) denote the den- 



*This article, except for Appendix B, is an improved version of 
Chapter IV of the author's Ph.D. dissertation, which has been sub- 
mitted in partial fulfillment of the requirements for the Degree of 
Doctor of Philosophy in the Department of Chemistry at Brown 
University. 



sity matrix for a similar free particle. 

The d-dimensional generalization of the Feynman-Kac 
formula is rather trivial. One just considers an inde- 
pendent Brownian bridge for each additional degree of 
freedom. To keep the notation simple, in this article 
we shall work exclusively with one-dimensional systems. 
However, the reader should notice that the main results 
of the paper remain true or have straightforward gener- 
alizations for systems of arbitrary dimensionality. 

In actual simulations, the Feynman-Kac formula is al- 
most always used in conjunction with Monte Carlo in- 
tegration methods^ and, for this purpose, one needs to 
construct rapidly convergent finite-dimensional approxi- 
mations to the stochastic integral described by Eq. (TJ). 
Ideally, such approximations should require knowledge 
of the potential only for the computation of the density 
matrix or the partition function of the physical system. 
This type of methods will be called direct methods. The 
main question we address in the present article concerns 
the rate of convergence of a class of discretization tech- 
niques as measured against the number of variables uti- 
lized for path parameterization. Throughout the paper, 
we assume that the potential V(x) is an infinitely diffcr- 
cntiable and bounded from below function. 

Until recently, the fastest direct method available (as 
order of convergence) has been the trapezoidal Trot- 
ter discrete path integral (DPI) methodi&£ The tech- 
nique is usually derived by means of the Lie- Trotter 
product formula and an appropriate short-time high- 
temperature approximation. The formal asymptotic con- 
vergence of the trapezoidal Trotter DPI method and 
of related DPI techniques was extensively studied by 
Suzuki^ and was found to be 0(l/ro 2 ). I shall com- 
ment more on this method in Section II. A. With the in- 
troduction of the random series implementation of the 
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Feynman-Kac formulaji . faster methods became avail- 
able. More precisely, two examples of direct path inte- 
gral techniques constructed around the Levy-Ciesielski 
and the Wiener-Fourier random series representations of 
the Brownian motion and pertaining to the general class 
of reweighted random series techniques were shown to 
have 0(l/n 3 ) asymptotic convergencei^^*^ In a recent 
Monte Carlo simulation^ the superior convergence of 
the reweighted methods proved to be crucial for the ac- 
curate determination of the potential, kinetic, and total 
energies of a highly quantum mechanical Lennard-Jones 
cluster made up of 22 molecules of hydrogen at a tem- 
perature of 6 K. 

In this article, I try to argue that in fact, for in- 
finitely differentiable potentials V(x), there might ex- 
ist direct short-time high-temperature approximations of 
arbitrary polynomial convergence order. The construc- 
tion of such approximations is based upon an "exper- 
imental" theorem on the pointwise convergence of Lie- 
Trotter product formulas. This theorem is presented in 
Section II. A, where it is used to derive a set of functional 
equations that short-time approximations must satisfy 
in order to have a given convergence order. Unlike stan- 
dard approaches based upon the construction of "effec- 
tive" potentialS(£*i2iiliii£iiLi£ the short-time approxima- 
tions we consider in the present article are based on 
carefully designed finite-dimensional approximations to 
the Brownian motion entering the Feynman-Kac formula. 
The potential itself is left unchanged. It is for this reason 
that the set of equations mentioned above do not depend 
upon the potential. The equations can be solved once 
for a given order and their (not unique) solutions can be 
tabled and used in actual computations for all potentials. 

The main mathematical problem that is left unsolved 
in this article is the existence of finite-dimensional ap- 
proximations to the Brownian motion that satisfy the 
functional equations for a given convergence order. To 
support the idea that such solutions exist, I explicitly 
construct two short-time approximations to the density 
matrix having convergence orders 3 and 4, respectively. 
A solution for the order 3 has been previously derivedjii 
but the one I construct in the present paper utilizes fewer 
path variables and fewer quadrature points. The solu- 
tion for the order 4 is derived as evidence that the gen- 
eral problem of constructing finite-dimensional approx- 
imations of arbitrary order is positively solvable. The 
fourth order method has numerical requirements simi- 
lar to the trapezoidal Trotter method (as ratio number 
of calls to the potential over number of path variables). 
The method has been recently utilized in the study of 
the heat capacity of the Nei3 cluster^ 

In Section V, I verify by numerical simulations the 
asymptotic convergence of the two short-time approx- 
imations discussed above. The definite agreement with 
the theoretical predictions is interpreted as proof that the 
theoretical development in the present article is mathe- 
matically sound. I conclude the paper by speculating 
that sequences of short-time approximations for increas- 



ing convergence orders (if they exist) may provide expo- 
nentially fast approximations for imaginary-time "prop- 
agated" wavefunctions, as measured against the number 
of path variables. I then analyze the implications of this 
hypothesis with respect to the solvability of the dynami- 
cal sign problem for real-time Feynman path integrals on 
a classical computer. 

In Appendix B, I derive the convergence constant for 
the celebrated trapezoidal Trotter path integral tech- 
nique. The convergence constant is verified by numer- 
ical simulations. The excellent agreement between the- 
ory and simulation is interpreted as further evidence that 
Theorem 2 is a valid mathematical statement (perhaps 
after further restrictions on its hypothesis). 



II. PRODUCT APPROXIMATIONS 

In the first part of this section, I review the classi- 
cal results of Suzuki concerning the order of convergence 
of a special family of short-time approximations. These 
results serve illustrate the main difficulties regarding the 
construction of short-time approximations having conver- 
gence orders higher than 2. I then state a theorem con- 
cerning the pointwise convergence of Lie- Trotter prod- 
uct formulas and discuss its implications with respect to 
the design of short-time approximations having superior 
convergence orders. In Section II. B, I introduce a spe- 
cial class of short-time approximations constructed by re- 
placing the Brownian motion appearing in the Feynman- 
Kac formula with appropriate finite-dimensional Gaus- 
sian processes. The functions utilized in the construc- 
tion of these finite-dimensional Gaussian processes will 
become the unknown variables for the systems of func- 
tional equations controlling the orders of convergence of 
the associated short-time approximations. These systems 
of functional equations are derived in Section III. 



A. A convergence theorem for product formulas 

One of the most fruitful approaches to constructing 
finite-dimensional approximations to the quantum me- 
chanical density matrix was given by Trotter £ It exploits 
the fact that {e~" H ;(3 > 0} is a semigroup of operators 
on L 2 (R), so that 



(2) 



or, in coordinate representation, 



dz{x\e~^ H \z){z\e-^ H \x'). (3) 



(In this work, the Hamiltonian, the kinetic operator, and 
the potential operator are denoted by the symbols H , K, 
and V, respectively.) The Trotter approximation theo- 
rem states that 



-f m = lim 



e -PK/n e -f)V/r, 



3 



in the sense of strong operator convergence. The quantity 



e -f3K/n e -l3V/n 



is called a short-time high-temperature approximation of 
the exact density matrix operator e~^ H / n . 

There has been a lot of research on the rate of conver- 
gence of the above approximation or of similar Trotter- 
like formulas. Of particular significance is Suzuki's 
workj& which treats the more general problem based on 
short-time approximations of the form 



-p(K+v) 



e -a /3V e -ax0V _ 



- hf)K e- aifw [l + 0{l3 v+1 )] 1 (4) 



where the sequences of non-negative real numbers 
ao, a±, . . . , ai and b%, b%, . . . , bi are palindromic and sum 
to 1. Following Suzuki, a short-time approximation 
f v {K, V;/3) is called of order v if 

e~K K+v > = f v (K,V;/3)[l + 0(P v+1 )}. 
In this casej£ 



-0{K+V) 



fu K, V; 



1 + 



v+l 



(5) 



[To be rigorous, I mention that Eq. (J5J has been proved 
for bounded operators A and B. The respective theorem 
states that the operator norm error of the final n-term 
Lie- Trotter product formula decays as fast as 1 jn v . How- 
ever, experience shows that the orders of convergence are 
correctly predicted even for the unbounded operators K 
and V . Moreover, if the non-existence theorem discussed 
below is true for bounded operators, then it is also true 
for the more general class of unbounded operators.] 

The more general splitting formula given by Eq. |4"| 
was considered by Suzuki in order to produce path inte- 
gral methods having faster asymptotic convergence. Un- 
fortunately, the following theorem of Suzuki (see Theo- 
rem 3 of Ref. 0) says that 

Theorem 1 (Suzuki nonexistence theorem) There 
are no finite-length splitting formulae ^ of order 3 or 
more such that the coefficients ao,b\,ai,... are all real 
and positive. 

The Suzuki nonexistence theorem limits the asymptotic 
order of convergence of this type of discrete path integral 
methods to 2, order of convergence that is attained for 
the following symmetric Trotter-Suzuki short-time ap- 
proximation 



-I3(K+V) 



■ ±0V e -0K e - 



^ y [l + 0(/3 3 )] (6) 



(or the one obtained by permuting V with K). 

The Suzuki nonexistence theorem serves to illustrate 
the difficulty of constructing path integral methods hav- 
ing asymptotic convergence better than 0(l/n 2 ). The 
idea of the Trotter theorem is commonly employed in 



the physical and chemical literature in order to generate 
faster integral methods starting with more general short- 
time approximations. The general strategy is as follows. 
Based upon a certain physical model, one constructs a 
short-time approximation po(x, x'; j3) of the true density 
matrix. Then, one corrects upon the short-time approxi- 
mation with the help of the Lie- Trotter product formula 



p n (x,x';(3) = [ dxi... [ dx n po (x,xi; - 
Jr Jr V n+1 

j. P 



Po x n ,x 



71 + 1 



(7) 



If the short-time approximation po(x, x'; f3) is "better" 
than the trapezoidal Trotter-Suzuki one, improved n-th 
order approximations to the exact density matrix may be 
obtained. The notion of "better" approximation may re- 
fer not only to the order of the short-time approximation 
but also to the overall quality of the approximation for 
finite n£ 

At this point, we remark that working with conver- 
gence theorems in operator norm is difficult and not par- 
ticularly helpful for actual developments of better short- 
time approximations. Indeed, the short-time approxima- 
tions are usually constructed in the configuration space as 
symmetric integral kernels po(x, x'; /3) and many proper- 
ties related to the norm operator topology are not readily 
available. Therefore, it is generally more convenient to 
use pointwise [in the space K 2 x [0, oo) of triplets (x, x'; /?)] 
convergence theorems of the type shown by the follow- 
ing theorem, which applies provided that po(x, x'; (3) is 
symmetric. 

Theorem 2 ("experimental") Assume that there ex- 
ists the linear (automatically Hermitian) operator T v i\), 
called a convergence operator, that associates to each in- 
finitely differ entiable and compactly supported function 
ipix) the square integrable function 



(Tvip)(x) = lim 



Then 



J R [pa(x, x';(3)- p(x, x'; (3))ip(x')dx' 



v+l 



(8) 



(9) 



lim (n + I)" [p n (x, x'; (3) - p{x, x'; 0)] = 
jf 1 (x \e-^ H T v e-^ 3H \ x') dd, 

where p n (x,x';/3) is defined by Eq. Q). 
Justification. Let T' v (x,x l '; 0) be defined such that 
p (x, x'; (3) = p(x, x';(3)+ P v+1 T' v (x, x>;0). 



Lie- Trotter composing the above relation n times and us- 
ing the semi-group property of the exact density matrix, 
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one argues that 



B. A general class of short-time approximations 



gv+1 r r 

p n (x, x'\ f3) = p(x, x'\ (3) + - 2_, dx ^ dx 

\^ ~* ) - n tJM. J M, 



j=0' 

xp \x,xi\ — — ) T„ \ xi,x 2 ; 



Xp \X2, X' 



n + 1 / \ n + 1 

,. (n-j)P\ , 



n + 1 



+ 0{l/n v+L ). 



In the limit n — > oo, one uses Eq. (|SJ) to cast the previous 
equation into 



lim {n + l) v [p n (x, x';{3) — p(x, x'; /?)] = lim 



Y^Jd XlP u Xl ,-^-\{T vP ) u 




In the formula above, the operator T v acts upon the den- 
sity matrix to the right through the first variable. Finally, 
one notices that in the same limit n — > oo, the Riemann 
sum transforms into an integral over the interval [0, 1], 
so that 

lim (n + 1)" [ Pn (x, x'; (3) - p(x, x'; 13)] = (3 U+1 
n — >oo 

x f d6 [ dx 1 p{x,x 1 ;9(3){T v p){x 1 ,x';{l-O)0) . 
Jo it 

Of course, Eq. © is nothing else than the above identity 
in Dirac's bra-ket notation. □ 
Observation. It is needless to say that the above justifi- 
cation is not a proof, nor is the hypothesis of the theorem 
completely stated. In keeping with the scope of the pa- 
per (and with the level of mathematical knowledge of the 
author), I only provide the basic reasons why the theo- 
rem must hold. The main effort of the present work is 
toward justifying the need for the theorems presented. 
The hope is that the mathematician will find the theo- 
rems interesting and worth investigating. However, all 
results deduced from this theorem, including the conver- 
gence constant for the trapezoidal Trotter path integral 
method (see Appendix B), are verified by numerical sim- 
ulations. 

Theorem [S] facilitates the construction of more accu- 
rate short-time approximations because it provides the 
exact convergence constant of the respective path inte- 
gral method in coordinate representation. In general, 
for a given order v, one would like to design short-time 
approximations po(x, x'; 0) that minimize (as absolute 
value) the convergence constant. In the ideal situation 
that the convergence constant is canceled, the order of 
convergence increases by one. In Section III, we shall 
use Theorem [5] to derive the set of equations that must 
be satisfied by the short-time approximations of a given 
order v. 



To make optimal use of Theorem [21 we need to de- 
vise systematic ways of constructing symmetric and pos- 
itive short-time approximations po(x,x';(3) for any or- 
der v. The positivity of the short-time approximation 
po(x,x';/3) is necessary in order to avoid the appear- 
ance of the sign problem in the Monte Carlo simulations. 
Development of such systematic ways has been previ- 
ously attempted by SuzukiSS as well as by Makri and 
Miller, 21 among othersii^ 2 ^ Unfortunately, all short- 
time approximations constructed so far involve deriva- 
tives of the potential V(x), derivatives that are either 
considered explicitly or introduced through the utiliza- 
tion of commutators involving the kinetic and potential 
operators. In fact, the higher the convergence order, 
the higher the order of the derivatives that are neces- 
sary. For this reason, except for the Takahashi-Imada 
approximation^ such approaches have enjoyed only lim- 
ited use. As discussed in the introduction, direct short- 
time high-temperature approximations based solely on 
the use of the potential function are more desirable. 

In this subsection, I present an alternative approach to 
constructing direct short-time approximations, approach 
that is related to the random series representation of the 
Brownian motion^ Evidence that will be presented in 
the subsequent sections supports the claim that the ap- 
proach is general enough to accommodate any arbitrary 
convergence order v. In this work, unless otherwise spec- 
ified, ao,a%,... denotes an infinite sequence of indepen- 
dent identically distributed (i.i.d.) standard normal vari- 
ables. The short-time approximations are constructed by 
replacing the Brownian motion in the Feynman-Kac for- 
mula with the finite dimensional Gaussian process 



y^QfcAfc(t 



(10) 



k=0 



The continuous and piecewise smooth functions 
{A k (u);0 < k < q} are required to satisfy the following 
relations: 



A o (0) = 0, A (l) = 1, and 

A fe (0) =A*(1) = 0, forl<fc<q. 



(11) 



The general expression of the short-time approximations 
we study in the present paper is 

p (x,x';(3) = pf p (x,x';/3) / dp,{a{) ■ ■ ■ / dp(a q ) 



x exp < — /3 J V 



k=l 



du), (12) 



where 



dp,(a k ) = (27r) 1/2 exp (-a 2 /2) da k 



and where 



x r (u) = x + (x' — x)Aq(u) 



5 



is a reference path connecting the points x and x' . 

A second condition we enforce on the system of func- 
tions {Afc(u); < k < q} is that 



Ao(«) + Ao(l-u) = 1 



(13) 



and that the finite dimensional Gaussian process 
Sfe=i a kA k (u) is invariant under the transformation u! 
1 That is, we require that 



B° u = ]T a k A k (u) i £ a fc A fe (l - u) = B?_ u . (14) 



fc=i 



fc=i 



The property described by Eq. (|14|l is analogous to the 
time symmetry of the standard Brownian bridge B®, 
which is the fact that {B®_ u , < u < 1} is also a Brown- 
ian bridge and is equal in distribution to {B®, < u < 1}. 
As a direct consequence of Eqs. l|T3|l and ifT^jl . the short- 
time approximation po(x, x'; f3) given by Eq. (|12fl is sym- 
metric under the permutation of the variables x and 
x 1 . This can be verified by performing the substitution 
v! = 1 — u in Eq. (|12fl . The time symmetry of the finite 
Gaussian process 53fe=i a fcAfc(u) can be enforced, for ex- 
ample, by restricting the functions {A k (u); 1 < k < q] to 
the class of symmetric and antisymmetric functions. 

In this general setting, given a fixed integer q, The- 
orem |21 suggests that the functions {A k (u);0 < k < q} 
should be chosen such that the order of convergence be 
maximized. We shall show in the next section that the 
system of functional equations controlling the order of 
convergence is independent of the potential V(x). This 
system of equations does not uniquely determine the 
functions {A k (u);0 < k < q}. For instance, it is a trivial 
matter to show that the short-time approximation given 
by Eq. 112fl is invariant under a linear orthogonal trans- 
formation of the functions {A k (u); 1 < k < q). 

A consequence of the constraint given by Eq. I|ll|) is 
the fact that the distributions of the end points B\ and 
B\ are identical and equal to that of the variable a§. 
In order to reproduce in a better way the properties of 
the Brownian motion, we may also require (but it is not 
necessary) that the pairs of Gaussian variables [B\,M\] 
and (Si, Mi) have equal joint distribution. Here, M\ and 
Mi are the so-called path centroids^ (first moments of 
the Brownian motion and its short-time approximation) 
and are defined by the equations 



Mi = / B u du 7 



Mi = I B u du and 
>o 



respectively. To find the class of short-time approxi- 
mations for which this condition is "built in," consider 
\q{u) — 1 and \i(u) = \/3(l — 2u), the first two nor- 
malized Legendre polynomials on the interval [0, 1]. Let 
{Afc(u)}fc>2 be a set of functions which together with the 
first two Legendre polynomials make up an orthonormal 
set on [0,1]. The Ito-Nisio theorem (see Theorem [fj] of 



Appendix A) says that 

oo 

B u = a u + ciiV3u(l - u) + a k A k (u), 

k=2 

where 

A fc (u) = / \ k (r)dT. 



Let us notice that if k > 2, then A^(l) = [by the 
orthogonality of (u) on 1] and 



K k (u)du = Afc(l) 



\ k (u)udu = 



[by the orthogonality of X k (u) on u). Therefore, Bi = do 
and 

1 V3 
Mi = -a Q + — ai 

depend solely on the variables aQ and ai. A little thought 
shows that we can build in the correct joint distribution 
of the end point and the path centroid by further re- 
stricting the class of functions {A k (u);0 < k < q} to 
those satisfying the constraints 



Ak(u)du = 1/2, forfc = 0, 
Jg 1 A k (u)du = x/3/6, for k = 1, and (15) 



jg 1 A k (u)du = 0, 



for 2<k<q. 



Until now, we have assumed that the path averages of 
the type 



V 



.(u) + ay^a k A k (u) 



k=l 



du 



are evaluated exactly. For practical applications, one also 
needs to devise a minimalist quadrature scheme specified 
by some points 0<ui < U2 < ... < u n < 1 and non- 
negative weights wi , W2 , • ■ ■ , w nq such that the discrete 
short-time approximation 

p (x,x';f3) = p f p(x,x';(3) / dp(ai) ■ ■ ■ / dp(a q ) 



x exp ■ 



i=i 



x r (u. t ) + Q-y^a fc A fc (uj) 
fe=i 



(16) 



has the desired convergence order. In the case of the 
discrete approximations, the set of quadrature points and 
weights as well as the values of the functions A k (u) at the 
quadrature points are fitting parameters. For the reason 
of ensuring time symmetry of the discrete formula, the 
quadrature scheme is required to be symmetric, i.e., the 
sequences 1/1,1*2 — tti, . . . , 1 — u nq and wi , W2, ■ ■ ■ , w Uq 
must be palindromic. 
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III. POWER SERIES EXPANSION FOR 
IMAGINARY-TIME PROPAGATED 
WAVEFUNCTIONS 

In this section, we shall derive the system of functional 
equations that must be satisfied by the functions Afc(«) 
appearing in Eq. (|1(J|> in order for the associated short- 
time approximation to have a convergence order v. To 
settle some terminology related to the utilization of the 
term "short-time," we interpret the parameter [3 as a 
time variable (physically, h(3 has dimension of time) so 
that the density matrix p(x,x';0) constitutes the time- 
dependent Green's function of a diffusion equation, or 
imaginary-time Schrodinger equation. As Theorem [21 il- 
lustrates, it is necessary to establish the power series ex- 
pansion of the imaginary-time propagated wavefunctions 
for the exact and the approximate propagators, respec- 
tively. I warn the reader that the power series derived 
in the present section are only a bookkeeping device for 
derivatives against f3 and are not required to converge to 
the actual imaginary-time propagated solutions. More- 
over, the potential V(x) and its derivatives are required 
to have finite Gaussian transforms. Actually, we require 
that 



1 



\j2ixa 



2 i^\v^(x + z)y 



< oo 



(17) 



for all x € M and a > 0, as well as for all integers k > 
and j > 1. This condition is necessary in order to en- 
sure that we recover the original potentials, derivatives, 
or products of such functions from their Gaussian trans- 
forms, in the limit that a — > (see Theorem 3 of Ref . 124}) . 



A. The exact propagator 

The power series expansion of the propagated wave- 
function 



is of utmost interest for the present development. With 
the help of the Feynman-Kac formula [according to 
Eq. (|A3|1 of Appendix A] and the Taylor power series 
expansion, one writes 



(x \e~^ H \ ip) =E \e-^o v(x+aB u )du^ x + aBi) 



E ■ 



oo 1 



j=0 ■ 



fe=0 



where 



{B u ) k du. 



(19) 



A second Taylor expansion leads to 



{x\e-P H \l>) = 

OO 



oo 1 



3=0 ■ 



k=0 



3=0 



We now expand the product in the preceding formula and 
collect the coefficients corresponding to the same power 

of p. Remembering that a — (fi 2 /mo)^ 2 fi 1 ^ 2 , one ar- 
gues that the powers of (3 are of the form where \i 
is a non-negative half integer, i.e. an element of the set 
N 2 = {n/2 : n € N}. For each /i£N 2 ,/i>0, define 



€ N 2m : ^2 kjk = 2/i [ . (20) 
fe=i J 



(x \e' pH | ip) = p(x,x';f3)i/j(x')dx' (18) A little thought shows that 



-PH 



i>) = E p 1 E (-iy hl ■■■■»■» (h 2 / mo ) 



31+33 + ^34+ ■■■ + (2^-2» 2m 



^eN 2 (j'i,...,j2 f .)eJ r f , 



^\x)V{x)^ [V^(x)] ... [V^- 2 ^(x)y 
ji!j 2 !...j 2/ ,!(2!)^(3!)^...[(2//-2)!]^ 



-E [(B^ (Mo)*" (Mip . . . (M 2AI _ 2 )^] , (21) 



with the convention that the term for [i = is ip{ x )- The 
fact that B u is a Gaussian distributed variable of mean 
zero implies that if j\ +js + 2j± + . . . + (2// — 2)j 2fJ , is odd, 
then 

E [(B^iMoy^A'hY 3 . . . (M 2 ^ 2 y*»] = 0, 



as can be verified by induction. Since ji + j 3 + 2j^ + 

. . . + (2/i - 2)j 2fl = 2fi- 2(j 2 H h j 2li ), one sees that 

h + h + 2j± + ... + (2/i - 2)j 2fJ/ is odd if and only if 
2/i is an odd integer. Thus, the sum in Eq. (|21J) can be 
restricted to the numbers fi € N 2 for which 2/i is even i.e., 
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the sum can be restricted to the set of natural numbers 



Therefore, 



p(x,x';P)ip(x')dx' = ^2/3" 53 (-l) j2+ - +h " (H 2 /m ) 
ip(M(x)V(xy* [vW(x)] h . . . [V^- 2 \x)\ ' 



M-(j2H hjap 



iiy 2 !...i2 |t !(2!)«(3!)«...[(2/i-2)!]^ 



E [(BiP(Mo)^(Mi)* . . . (M 2M - 2 ) i3 "] 



r 



(22) 



Observation. The power series expansion for the 
imaginary-time propagated wavefunction can also be de- 
rived by expanding the operator e~@ H in a power series. 
One obtains 

(x\e-? H \^) = Y,-{x\{-mf\Tp) 



oo 

5>4 



M=0 



ft 2 ri 2 

2mo c?x 2 



VW- (23) 



Of course, the terms of the two series given by Eqs. 11' I'D 
and \\'2'M are equal. However, as we shall see in the 
following subsection, Eq. (|22|) applies in an almost un- 
changed form for all short-time approximations defined 



by Eqs. I|12|l and (|16|) . In contrast, there might be no 
formal analogue of Eq. ll'.'ill for such short-time approxi- 
mations. 



B. The approximate propagator and the identities 
controlling its order of convergence 

The only property used for the derivation of the power 
series expansion of the exact propagator was the fact that 
the Brownian motion is a Gaussian process. Since the 
approximation to the Brownian motion given by Eq. i|10|) 
is also a Gaussian process, Eq. H22JI remains true for the 
approximate propagator, too. Therefore, 



oo 

/ p (x, x'; p)ip(x')dx' =5^/3" ^ (_!)*+...+*, (H 2 /m ) 
R a*=o (j 1 ,...,i 2fl )eJ (1 

^(x)v( x yi [v {1 Hx)Y 3 ■ ■ ■ [v {2 »- 2) (x)Y 21 



M-O'aH 



./;!.„!... .y 2 „!;2!).'M:i!)' ... (2// 2)!]^ 



(B 1 )^(M ) J2 (M 1 



(24) 



where 



M k = 



du. 



(25) 



If the discrete short-time approximation given by 
Eq. I|16fl is employed, then Eq. I|24(l remains true pro- 
vided that M k is redefined to be 



"5 j 



(26) 



if and only if 



= E 



(B^ (M y 2 (Mx)^ . . . (M 2m _ 2 )^ (27) 



for all 2/j-tuples of non-negative integers (ji, j 2 , • • ■ , J2^) 
such that 

53 k 3k = 2p 

fc=i 

and 1 < (i < v. 



Theorem [5] immediately implies the following state- 
ment. 



The general problem that one would like to solve us- 
ing the theory developed so far is the following. Given 
a convergence order is there a finite system of func- 
Theorem 3 A short-time approximation of the types tions A k (u) such that the corresponding short-time ap- 
given by Eq. \1U\) or Eq. \lb\) has convergence order v proximation has order vl If the answer is yes, what is 
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the minimal number q of functions A/-(it) necessary to 
achieve the respective convergence order? Then, what is 
the minimal number of quadrature points such that a dis- 
crete short-time approximation has convergence order vl 
The relevance of the questions asked in the current para- 
graph will be further clarified in Section VI, where we 
analyze the problem of minimizing the statistical noise 
for real-time propagators. 

IV. EXAMPLES OF SHORT-TIME 
APPROXIMATIONS HAVING CONVERGENCE 
ORDER 3 OR 4 

In this section, I try to present evidence in support of 
the idea that the system of equations appearing in The- 
orem for a given order v is always satisfied by some 
finite system of functions A^(m). I do this by comput- 
ing explicit numerical solutions for the convergence or- 
ders 3 and 4. As apparent from Table U the number of 
equations that need to be verified for a given order v in- 
creases rapidly with v. In fact, the number of elements 
of J M is the number of distinct partitions of 2fi. With 
the help of the Hardy-Ramanujan asymptotic formula, 27 
one deduces that the number of equations that need to 
be verified for a given order v behaves asymptotically as 

Therefore, the "by hand" approach utilized in the present 
section is bound to fail even for slightly larger conver- 
gence orders. By use of computers, one may hope to ob- 
tain solutions for moderately large convergence orders. 
However, I believe future work on the problem may re- 
veal better strategies for the computation of short-time 
approximations of high convergence orders. 

The two short-time approximations constructed in 
the present section are called reweighted short-time 
approximations. — The defining features are the equality 
Ao(it) = u and the fact that the functions {Ak(u);l < 
k < q} appearing in Eq. ((lTj|l are required to satisfy the 
constraint 

<? 

J2h(u) 2 = u(l-u). (28) 
fe=i 



The last equation stems from the condition that the 
Gaussian variables B u and B u have equal variances for 
each u £ [0,1] (equal weights). As we shall see, if this 
constraint is imposed, most of the functional equations 
for convergence orders 3 and 4 are automatically satis- 
fied. However, the number of remaining equations still 
scales exponentially and, for higher convergence orders, 
the constraint given by Eq. I|28|l may actually become a 
nuisance. 

One additional feature of the reweighted short-time ap- 
proximations stems from the relation Aq(u) — u and fa- 
cilitates the numerical implementation of the associated 
Lie- Trotter product formula given by Eq. (J7J). The fol- 
lowing generalization of a result of Predescu and Doll (see 
Theorem 2 of Ref. I23T) is straightforward to prove. 

Assume n is of the form n = 2 fc — 1 and let {a;,-; 1 < I < 
k, 1 < j < 2 1 ' 1 } and 1 < I < q, 1 < j < 2 k } be two 
independent sets of i.i.d. standard normal variables. Let 
{Fij(u);l > 1,1 < j < 2 l ~ x } be the system of Schauder 
functions^ on the interval [0, 1]. The Schauder functions 
can be generated by translations and dilatations as fol- 
lows. Let F x>l (u) : R -> K be defined by 

f u, ue [0,1/2], 
F hl {u) = < 1-u, u 6(1/2,1], (29) 
I 0, elsewhere. 

Then, 

F ld (u) = 2-«- 1 )/ 2 F 14 (2'~ 1 U - j + 1), (30) 

for all Z > 1 and 1 < j < 2 l ~ 1 . Extend the functions 
{A;(m);1 < I < q} outside the interval [0,1] by setting 
them to zero [the same way the first Schauder function 
Fi,i(u) was extended to the whole real axis in Eq. 12911 ] 
and define 

Gi,j{u) = 2- k / 2 A l (2 k u- J + l), (31) 

for 1 < / < q and 1 < j < 2 k . 

In these conditions, the following theorem holds. 
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Theorem 4 With the convention that a ; 2 '- 1 +i — an d \ 2 fc +i = for all I G 1, k, we have 

k 2'- 1 



p n {x,x';f3) 
p fp (x,x';/3) 



dais... I da k>2 k-i (2tt) " /2 exp[-^^^ 

i=i j=i 



1 



a h 



x / dfri 



x exp < —j3 / V 



db, 2k ( 27 r)-(" +1 ^ 2 exp 

\ i=i i=i 

a /,[2'- 1 «] + l Fl,[2'- 1 u] + l( u ) 



(32) 



2=1 



"<7 y^&j,[2»ttl+l G , ; ! [ 2 fc u] + l( M ) 
I 



/ = 1 



where [2 l 1 u] and [2 fe w] are the integer parts of2 l 1 u and 
2 k u, respectively. 

The reader can easily verify that Eq. IMl'l) is a so-called 
reweighted Levy-Ciesielski path integral technique, as de- 
fined in Ref . [Ill It has been argued^ that this represen- 
tation is more advantageous than the direct expression of 
p n {x, x '; (3) that is obtained from the Lie- Trotter product 
formula, for practical implementations. The expression 
obtained by Lie- Trotter composing the discrete version of 
Po{x, x'; (3) given be Eq. I|16[l can also be put in the form 
of Eq. ll.j2H . However, the one-dimensional integral at ex- 
ponent is replaced by a quadrature sum. The quadrature 
scheme is specified by the n q 2 k (not necessarily different) 
quadrature points 

u'u = 2- k { Ul +j - 1), 1 < i < n q , 1 < j < 2 k (33) 

and the corresponding weights 



w' id = 2 k Wl 



(34) 



The new quadrature points u[ ^ are obtained by trans- 
lations and dilatations (more precisely, contractions) of 
the original quadrature points U{. 



A. Reweighted short-time approximation having 
convergence order 3 

The equations that the functions A/~(u) must satisfy 
in order to generate a reweighted short-time approxima- 
tion of order 3 are those of the type shown by Eq. I|27|) 
for the indexes {31,32, ■ ■ ■ , j2/x) presented in Table [IJ with 
p, = 1, 2, and 3. For a better understanding, we men- 
tion that in Table ^ we only present the non-zero com- 
ponents of a given index (ji, 32, ■ ■ ■ , J2^)- There are a 
total of 2 + 5+11 = 18 equations that should be veri- 
fied. However, given the special form of the reweighted 
finite-dimensional approximation to the Brownian mo- 
tion, most of these equations are automatically satisfied. 



du 



TABLE I: Indexes of the equations that need to be verified 
for various values of p. Shown are the non-zero components 
of these indexes. 





= 1 




32 = 1 






js = 1 
















37 = 1, ji = 1 
















je = 1, 32 = 1 










J4 = l 






je = 1, ji = 2 










h = 1, ji = 1 






js = 1, 33 = 1 




/' 


= 2 




32 = 2 




js 


= L 32 = 1, ji = 


1 








32 = 1, ji = 2 






js = 1, ji = 3 










ji = 4 






ji =2 














3* 


= 1, j3 = 1, jl = 


1 














ji = 1, 32 = 2 










je = 1 


p = 4 


Ji 


= 1, j2 = 1, jl = 


2 








35 = 1, ji = 1 






ji = 1, ji = 4 










ji = 1, j2 = 1 






33 = 2, j 2 = 1 










34 = 1, ji = 2 






j 3 = 2, ji = 2 










33 = 1 




js 


= 1, j 2 = 2, ji = 


1 


/' 


= 3 


33 


= 1, .72 = 1, ji = 1 




33 


= L 32 = 1, ji = 


3 








33 = 1, ji = 3 






j3 = 1, jl = 5 










32 = 3 






j 2 =4 










ji = 2, ji = 2 






j 2 = 3, ji = 2 










32 = L ji = 4 






32 = 2, ji = 4 










ji = 6 






j2 = 1, jl = 6 
















ji =8 





As such, the equations for which the only non-zero com- 
ponents are j\ and j 2 are verified by all reweighted short- 
time approximations. The discrete versions satisfy the 
respective equations provided that 



E 



IV i 



1. 



One actually checks that all equations for p = 2 as well 
as all equations for p — 3, except for the one specified 
by 32 = 2, are automatically satisfied. The discrete ver- 
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sion verifies these equations provided that the quadrature 
scheme is capable of integrating exactly all polynomials 
1, u, and u 2 . For example, let us consider the equation 
specified by j§ = 1 . We have 



E 



E w *(^) 4 



,i=l 



E "•< E f(^) 4 l = 3 E 



i=l 



By Eq. (|27[) as specialized for j 6 = 1, the above value 
should equal [see Eq. I|A2|1 of Appendix A] 
fi 



(B u fdu 



f E [(B u ) 4 ] du = 3 / 
o Jo 



u du. 



This shows that the quadrature technique must integrate 
exactly the polynomial u 2 . 

We now turn our attention to the remaining equation 
defined by j'3 = 2. One computes 



(35) 







2 




2 q 




E 


/ B u du 




/ udu 


+E 


/ A k (u)du 




Jo 




Jo 


fc=i 


Jo 



which should equal 







2 




2 




E 


/ B„ du 
Jo 




/ udu 
Jo 


+ 3 


[L 



l 2 

u(l — u)du 

To compute the expected value of the square of the first 
moment of the Brownian motion, write the Brownian mo- 
tion as a random series constructed via the Ito-Nisio the- 
orem from the Legendre orthogonal polynomials on the 
interval [0,1]. Then, as discussed in the preceding sec- 
tion, 

/ B u du = ciq udu + VSai / u(l — u)du 
Jo Jo Jo 

and Eq. (|36|l follows. From Eqs. I|35|l and l|36|l . one easily 
obtains the identity 

, 2 

Kk{u)du 



1 v l-l 



E 

fc=i 



1 

12' 



A similar relation can be deduced for the discrete ver- 
sion but with the integrals replaced by the corresponding 
quadrature sums. 

We can summarize the findings of the present subsec- 
tion into the following proposition. 

Proposition 1 A reweighted short-time approximation 
has order 3 if and only if 

q r „1 n2 



E 

fc=l 



Kk{u)du 



Uo 



1 

12' 



(37) 



A discrete reweighted short-time approximation has or- 
der 3 provided that the associated quadrature scheme in- 
tegrates exactly all polynomials of degree at most 2 and 
provided that 

9 [«! "I 2 1 

E E m ^ fe ( ui ) 

fc=l Li=l 



12' 



(38) 



We conclude the present subsection by constructing 
a minimalist reweighted short-time approximation hav- 
ing convergence order 3. Because of the identity 
the minimal number q of functions A k (it) capable of sat- 
isfying Eq. (|37|l is 2. Indeed, if q = 1, then Ai (it) = 
[u(l - u)} 1 / 2 and 



1 ~ I 2 
Ai(u)du 



11 



tt 2 /64^ 1/12. 



We now try a set of two functions of the form 



Ai(u) = yju(l — u) cos[a(u — 0.5)], 
^{u) = \/u(l — u) sin[a(u — 0.5)]. 



(39) 



The functions Ai (it) and A2(u) are orthogonal because 
the first is symmetric under the transformation u' = 1—u, 
whereas the second is antisymmetric. The constant a 
is then determined by Eq. (|37|l and has been evaluated 
with the help of the Levenberg-Marquardt algorithm, as 
implemented in Mathcad^ The solution has the approx- 
imate value 



3.056620471. 



(40) 



To design a minimalist discrete short-time approxi- 
mation of order 3, we consider an arbitrary symmetric 
quadrature rule on the interval [0, 1] that integrates ex- 
actly all polynomials of degree less or equal to 2. Then, 
we find the value of a that satisfies Eq. (|38|l for the cho- 
sen quadrature technique. It is not difficult to argue that 
the minimal number of quadrature points in the open in- 
terval (0, 0.5) must be 1. The reason is that the values of 
the functions Ai(u) and A2(u) at the points u = and 
u = 0.5 do not depend upon the parameter a. Thus, 
Eq. I|38|) cannot be satisfied if there are no quadrature 
points located inside the open interval (0, 0.5). 

The quadrature rule is taken to be the 2-point Gauss- 
Legendre rule on the interval [0, 1], quadrature rule that 
integrates exactly all polynomials of degree less or equal 
to 3. The appropriate value for the parameter a is then 
determined from Eq. I|38|) and is found to be 



2.720699046. 



(41) 



The quadrature scheme is given in Table ITU for ease of 
reference. 



TABLE II: Quadrature points and weights for the minimalist 
discrete short-time approximation of order 3. The points and 
weights are those for the 2-point Gauss-Legendre rule on the 
interval [0,1]. 



i 


1 


2 


Ui 


0.211324865 


0.788675135 


Wi 


0.500000000 


0.500000000 



As shown by Eq. (|32|l . the number of path variables 
entering the expression of p n {x, x'; (3) is (q+l)n+q = 3n+ 
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2, whereas the number of quadrature points [see Eq. Q33[l] 
is n q (n + 1) = 2n + 2. Thus, for large enough n, the ratio 
(2n + 2) /{in + 2) approaches 2/3, value that is smaller 
than the one for the trapezoidal Trotter discrete path 
integral method. Therefore, the method described in the 
present paragraph has fewer numerical requirements than 
the trapezoidal Trotter discrete path integral method for 
equal numbers of path variables, yet it achieves cubic 
convergence for smooth enough potentials. 



B. Reweighted short-time approximation having 
convergence order 4 



Because the number of equations to be verified in- 
creases significantly for the reweighted short-time ap- 
proximations of order 4, we choose to approximate the 
Brownian motion by the finite dimensional process 



B u = a u + ai\fZu(l - u) + ^ a k A k (u), 



(42) 



k=2 



where the functions A k (u) satisfy the equations 



/ 

Jo 



A k (u)du = 0, for 2 < fc < q. 



As discussed in Section II.B, in this case the vari- 
ables (Bo, Mo, Mi) have the same joint distribution as 
(Bo, Mo, Mi) (notice that Mo and Mq are equal con- 
stants). This remains true of the discrete reweighted 
short-time approximations provided that the quadrature 
scheme integrates exactly the polynomials of degree at 
most 2 as well as the functions A k (u), for 2 < k < q. 

Using the special form of Eq. I|42|) , it is not difficult to 
verify that all the equations in Tableware automatically 
satisfied with the exception of the one specified by ji = 2. 
This remains true of the discrete versions provided that 
the quadrature scheme integrates exactly all polynomials 
of degree at most 3 as well as the functions A k (u), for 
2 < k < q. For the sake of an example, let us consider 
the equation specified by js = l,j3 = 1, which is the 
most difficult to verify. I leave it for the reader to argue 
that in general 



E E 

("il *2"i2 ^3 "14 

\ii,i2,i3M 
= E ( M M,M + M id,i,j + M i,j,3,i) ■ 



(43) 



Using Eq. (|43(l . one computes 

E^ ^ udu J Q " 3 E J Q M u ) du 



x / Ai(u)Aj(uYdu 
'0 



i,j=0 
1 v 1-1 



= 3E 



i=0 



Ai(u)du 



x / Ai(u)udu 
Jo 

x / Ai(u)udu 
Jo 



1 1 

2 + 8 

1 1 

2 + 8' 





<? r d 



3E 



i=2 



Ai(u)du 







where we used the equality 

The above equation remains true of the discrete versions, 
too. For the full Brownian motion, one computes via 
the random series representation based on the Legendre 
orthogonal polynomials on the interval [0, 1] 



and the fact that the equation j'5 = 1, j'3 = = 1 is satisfied 
follows. 

We now turn our attention to the equation specified 
by j?4 = 2. One computes 



1 x 2 
2, 



Btdu 



E 



E I ciidj 



du 



where 



-4,3 — I Ai(u)Aj(u)du. 
10 



Using Eq. one deduces 



E 



^° ' i,i=0 \i=0 



At this moment it is useful to remember that Ao(u) = u 
and A\(u) = \/3u(l — u). Moreover, notice that Eq. I|28() 
implies 

1 ,1 1 

y^^Cj.i = / [u 2 + u(l - u)]du = -. 

i—rt ^0 



Therefore, 
E 



1,3 



I 



I 2 

Bldu 



^E 



1 - I 2 

Ai(u)Aj(u)du 
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For the full Brownian motion, one computes via the 
Wiener-Fourier series 



E 



1 x 2 

2, 



Bidu 



fc=i 



- + 4V 
9 ^ 

fc=i 



1 ' 2 sin(/c7TM) 

M-l/ — ; du 



2 sin(/c7ru) 



1 

4 



OO 



7T 4 ^ fc 4 
k=l 



+ - 



1 

E p 



Then, the equality 
i- 1 



E 



k=l 



E 



A' 4 



1 

4 



1 N 2 

Bidu 



implies 



E 

i,3'=0 



Aj(tt)Aj-(u)du 



(44) 



With the one-dimensional integrals replaced by appropri- 
ate quadrature sums, Eq. (|44fl must also be satisfied by all 
discrete short-time approximations of order 4. Remem- 
ber that the quadrature scheme is assumed to integrate 
exactly all the polynomials of degree at most 3 and all 
the functions Ak(u) for 2 < k < q. 

In the remainder of this subsection, we construct an 
example of reweighted short-time approximation of order 
4. Clearly, we cannot set q = 2 in Eq. (|42|l because then 

A 2 (u) ^ {u(l - u)[l - 3u(l - u)}} 1/2 , 
as follows from Eq. (|28|l . and consequently, 



A 2 (u)du ^ 0. 

Thus, we set q — 3 and look for functions of the form 

!A 2 (u) = r(u) cos[ai(u — 0.5) + a%{u — 0.5) 3 ], , 
A 3 (u) = r(u) sin[ai(u - 0.5) + a 2 (u - 0.5) 3 ], 

where 



1/2 



r(u) = {u(l - u)[l - 3u(l - u)}} 



The functions A-2(u) and A3 (it) are orthogonal because 
the first is symmetric under the transformation v! = 1 — 
u, whereas the second is antisymmetric. The integral 
over [0, 1] of the function hz(u) is zero by antisymmetry. 
Then, the constants aj and a 2 are determined from the 
system of equations 



A 2 (u)du = 0, 




3 

E 

i,j=0 



1 _ _ -|2 

Ai(u)Aj(u)du 



(46) 



The values of the constants a% and a 2 have been deter- 
mined numerically to be 

ai w 5.768064999 and a 2 « 13.49214669. (47) 

Let us now design a minimalist discrete short-time ap- 
proximation of order 4. Given an arbitrary symmetric 
quadrature technique that integrates exactly all polyno- 
mials of degree less or equal to 3, we determine new values 
for a± and a 2 from the system of equations 



1) 



2) 



1=1 
3 

E 

i,j=0 



wiA 2 (ui) = 0, 



1=1 



wiAi(ui)Aj(ui) 



(48) 



Because there are two equations, it is easy to argue that 
the number of quadrature points lying in the open in- 
terval (0, 0.5) must be at least two. Consistent with this 
observation, the quadrature technique is chosen to be the 
4-point Gauss-Legendre technique on the interval [0,1]. 
This quadrature technique integrates exactly all the poly- 
nomials of degree at most 7. The new values for the 
parameters a.\ and a 2 are then determined by solving 
the system of equations given by Eq. I|48|l for the cho- 
sen quadrature scheme. The solution of the system of 
equations is given by 

ai re 6.379716466 and a 2 « 8.160188248. (49) 

The quadrature weights and points are presented in Ta- 
ble IIIII for ease of reference. 



TABLE III: Quadrature points and weights for the minimalist 
discrete short-time approximation of order 4. The quadrature 
points and weights are those for the 4-point Gauss-Legendre 
technique on the interval [0, 1]. 



i 


1 


2 


3 


4 


Ui 


0.069431844 


0.330009478 


0.669990522 


0.930568156 


Wi 


0.173927423 


0.326072577 


0.326072577 


0.173927423 



As shown by Eq. (|32|l , the number of path variables en- 
tering the expression of p n {x, x'; 0) is (q+l)n+q = 4n+3, 
whereas the number of quadrature points [see Eq. 1)33(1 ] 
is n q (n + 1) = An + 4. Thus, for large enough n, the ratio 
(An + 4)/ (An + 3) approaches 1, value that equals the one 
for the trapezoidal Trotter discrete path integral method. 
Therefore, the fourth order method has the same numer- 
ical requirements as the trapezoidal Trotter discrete path 
integral method for equal numbers of path variables, yet 
it achieves quartic convergence for smooth enough poten- 
tials. 



V. NUMERICAL VERIFICATION OF THE 
ASYMPTOTIC ORDERS OF CONVERGENCE 

One of the main advantages of the Lie- Trotter prod- 
uct formula consists of the fact that, for low dimensional 
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systems, the evaluation of the density matrix and related 
properties can be performed accurately by means of the 
numerical matrix multiplication (NMM) methodi^SiSi We 
shall use the NMM method to compute n-th order ap- 
proximations to the partition function of the type 



for one-dimensional systems. We follow closely the simu- 
lation strategy employed in Ref . ^2 for a similar numeri- 
cal study of asymptotic orders of convergence. The sym- 
bol (v) to the exponent serves to differentiate between 
short-time approximations of different orders v. 

The main steps of the NMM algorithm are as follows. 
First, one restricts the system to an interval [a, b] and 
considers a division of the interval of the type 

Xi = a + i(b - a) /M, < i < M. 

Next, one computes and stores the symmetric square ma- 
trix of entries 



a _ h ~ a M ( .I 3 



M 



n+l 



, 0<i,j<M. 



The value of the partition function can then be recovered 
as 

Z^([3)=t r (A n+1 ). 

By computer experimentation, the interval [a, b] and the 
size M of the division are chosen such that the compu- 
tation of the partition function is performed with the 
required accuracy. A fast computation of the powers 
of the matrix A can be achieved by exploiting the rule 
^rn+n _ (^jm^n _ p Qr more details, the reader is referred 
to the cited literature. 

The Gaussian integrals appearing in the expression of 
the discrete reweighted short-time approximation 

p$\x,x';/3) = p fp (x, x';(3) / d^( a i) ' • • / dfi(a q ) 



x exp 



E 



WiV 



E 

k=l 



a k k k (ui) 



can be evaluated by means of the Gauss-Hermite quadra- 
ture technique^! for small enough q (in our case, q is 2 
for the approximation of order 3 and 3 for the approxi- 
mation of order 4, respectively). For the purpose of es- 
tablishing the asymptotic convergence of the partition 
functions, it was found that a number of 10 quadrature 
points for each dimension is sufficient for both short-time 
approximations studied in the present section. This is so 
because the errors due to the Gauss-Hermite quadrature 
approximation quickly vanish as 13/ {n + 1) — > 0. 

Once the partition functions are evaluated, we com- 
pute the quantities 



R 



W) 

2m+ 



(50) 




o^a (3) ,-a (3) 



— 1 1 f~ 



10 20 30 40 50 
m 

FIG. 1: The convergence orders of the two discrete short- 
time approximations for the quartic potential. The plotting 
symbols are shown only for every tenth data point actually 
computed. 



and 



i^} = m 2 In 



As demonstrated in Ref. Hfl the slope of a™ as a function 
of to converges to the convergence order. We want to 
verify whether or not this convergence order is v. The 
exact partition function Z(J3) necessary in Eq. i|50|) is 
evaluated either by variational methods or by employing 
a large m. 

The first example studied is the quartic potential 
V(x) — x 4 /2. The following values of the physical 
constants (in atomic units) have been utilized: h = 1, 
777 = 1, and P = 10. The second example studied con- 
sists of a particle trapped on a line between two atoms 
separated by a distance Lm^ The particle is assumed to 
interact with the fixed atoms through pairwise Lennard- 
Jones potentials. The resulting cage is described by the 
potential 



V(x) = 4e 



(- 



L 



L 



if < x < L, and V(x) = +oo, otherwise. The parame- 
ters of the system are chosen to be those for the He atom. 

o 

We set TO = 4 amu, e/k B = 10.22 K, a = 2.556 A, and 

O 

L = 7.153 A. At T = 5.11 K, which is the temperature 
utilized in the present computations, the system is prac- 
tically in its ground state. For more details regarding 
the pres ent simulations, the reader is advised to consult 
Ref-IH 

As Figs. H an d El show, the orders of convergence pre- 
dicted in the preceding section are well verified. I inter- 
pret these results as proof that the mathematical analysis 
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FIG. 2: As in Fig. 0for the He cage problem. 



performed in the present paper is sound. The He cage 
problem is interesting because the Lennard- Jones poten- 
tial lies outside the class of potentials for which the the- 
ory was developed. As explained in Ref. Il2t the density 
matrix of the Lennard- Jones potential has an exponen- 
tial decay near singularities and therefore, the behavior 
of the potential near singularities is not important as far 
as the polynomial convergence of imaginary-time path 
integral methods is concerned. 



VI. CONCLUSIONS 

In this article, I have considered the problem of con- 
structing direct short-time approximations to the density 
matrix of a physical system of arbitrary convergence or- 
ders. I have shown that the problem can be reduced 
to the construction of finite-dimensional approximations 
to the Brownian motion that satisfy a certain system 
of functional equations. Using the developed theory, I 
have constructed two examples of reweighted short-time 
approximations having convergence orders 3 and 4, re- 
spectively. The predicted orders of convergence have 
been verified by numerical simulations. In Appendix B, I 
have derived the convergence constant for the trapezoidal 
Trotter path integral method. The predicted convergence 
constant has also been verified by numerical simulations. 

For imaginary-time path integral simulations, the 
reader may object that the use of a path integral tech- 
nique having faster asymptotic convergence is not a sig- 
nificant algorithmic improvement because the final com- 
putational effort is eventually controlled by the rate of 
convergence of the Monte Carlo integration method. The 
computational effort, as measured against the number of 
calls to the potential V(x), can be evaluated as follows. 
To attain a given absolute error e, one must utilize a 
number of 



path variables (here, const is some proportionality con- 
stant). The cost to evaluate the average potential for a 
given path is equal to the number of quadrature points, 
which, in turn, is proportional to the number of path 
variables [here, we do not take into account the cost for 
the computation of the paths, which scales as nlog 2 (?i), 
but which is usually negligible for the values of n com- 
monly employed in practice]. Thus, the cost for a single 
path evaluation is const/e 1 ^ . This cost is to be multi- 
plied by the number of Monte Carlo steps, which is given 
by the formula 

Nomc = const /e , 

assuming that the variance of the Monte Carlo method 
does not depend upon the number of path variables. 
Thus, the total cost, defined as the number of calls to 
the potential necessary to attain a given error, is 



Cost = const .£-(2+!/^ 



(51) 



where v is the convergence order of the direct path in- 
tegral method. Eq. I|51|l shows that we cannot beat the 
slow convergence of the Monte Carlo integration scheme 
by increasing the order of convergence of the path in- 
tegral technique. The total cost changes from e~ 2 - 5 to 
£ -2.25 oniy^ as we switch from the trapezoidal Trotter to 
the fourth order method designed in the present article. 

However, the methods designed in the present pa- 
per are still useful because the improvement, even if 
marginal, comes "free of any charge." Indeed, as shown 
in Section IV. B, the ratio number of quadrature points 
over number of path variables is 1 (for n large enough) 
for both the trapezoidal Trotter and the discrete fourth 
order method introduced in the present article. There- 
fore, there is no loss of efficiency in employing the discrete 
fourth order method even for those potentials for which 
the optimal convergence order is not attained. Because 
no additional cost is incurred even in the most disadvan- 
tageous situations, the discrete fourth order short-time 
approximation is a natural replacement for the trape- 
zoidal Trotter short-time approximation in all path inte- 
gral simulations. 

At a more general level, the present development may 
be relevant for the problem of performing real-time path 
integral simulations^ In this case, the asymptotic rate 
of convergence is crucial because the noise in the Monte 
Carlo simulation not only that does depend upon the 
number of path variables, but actually increases exponen- 
tially fast with the number of path variables. This is the 
statement of the well-known dynamical sign problem^ 

Let us assume that for a given convergence order v, 
there is a finite system of functions {Afc(u); < k < q v } 
that generates the short-time approximation of order v 



p$>(x,x';f3)=pf p (x,x';p) I dp, cr (a 1 ) ■ ■ ■ I da a {a q ) 



n = const /e 1 /" 



x exp < — J V 



•(u) + a k K k (u) 



k=l 



du 



(52) 
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where 



dfj, a (a,k) = (27n7 2 



-1/2 



exp [-a 2 k /(2a 2 )] da k . 



Notice that in Eq. (|52J) we have performed a substitu- 
tion of variables a' k = aat so that the dependence of the 
spread of the paths with (3 is no longer buried in the 
potential [remember, a = (fi^fl/mo) 1 / 2 ]. In principle, 
this transformation should allow us to extend the above 
formulas to complex-valued (3. We ask the question of 
whether or not it is more optimal to give up the use of 
the Lie- Trotter product formula altogether and instead 
consider the sequence of approximations 



p^\x,x';/3) -> p(x,x';/3) 



oo. 



(53) 



If with appropriate restrictions on V{x) and ip( x ) the se- 
ries appearing in Eq. j^H is analytic in (3, it is straight- 
forward to see that 



an alternative formulation of the Feynman-Kac formula 
and the random series construction of the Brownian mo- 
tion are presented. For further information, the reader is 
advised to consult the cited mathematical literaturei^2& 
Chapters I and II of Ref. [35J also contain an in-depth in- 
troduction to Brownian motion and its relation to the 
Feynman-Kac formula. 

A standard Brownian motion is defined as a stochastic 
process {B u ,u > 0} that satisfies the following condi- 
tions: 

(a) Given < uq < u\ < ... < u n an arbi- 
trary finite sequence of increasing times, the initial 
position B UQ and the position increments B Ul — 
B UQ , .... B Un — B Un l are independent. 

(b) If s, u > and [a, b) C K is some arbitrary interval, 
then 



p^\x,x'\p)i){x')dx' -> J p(x,x';f3)i/j(x')dx' (54) P {B u+S - B u e [a,b]) = exp dx. 



exponentially fast as measured against v. 

It is then apparent that a favorable scaling of q v with 
u, as for instance a polynomial scaling, may strongly 
alleviate the dynamical sign problem. As the Hardy- 
Ramanujan formula shows, the number of equations that 
must by satisfied by the system of functions {Afc(w); 1 < 
k < Qu} increases with v faster than any polynomial. 
However, this does not necessarily imply that q v increases 
with v at the same rate. In the examples constructed in 
Section IV, we have been able to accommodate the 18 
equations for order 3 with only two functions, whereas 
the 40 equations for order 4 were accommodated with 
three functions. In both cases, the actual number of func- 
tions was much lower than the number of equations. I 
hope this short analysis justifies my belief that future 
research on the subject is worth the time of investiga- 
tion and may lead to significant progress in the area of 
real-time path integral simulations. 
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APPENDIX A: SOME MATHEMATICAL FACTS 
ABOUT THE BROWNIAN MOTION AND THE 
FEYNMAN-KAC FORMULA 

In this appendix, I review the definition and some of 
the basic properties of the Brownian motion. In addition, 



(c) With probability one, the Brownian motion sam- 
pling paths B u are continuous. 

The existence of a stochastic process satisfying the above 
conditions has been first proved by Wiener— in 1923. 

If Bq = with probability one, then the Brownian 
motion is said to start at zero. In the present work, B u 
always denotes a standard Brownian motion starting at 
zero. The conditions (a) and (b) above are sufficient to 
demonstrate that the Brownian motion starting at zero 
is a Gaussian process with joint finite distributions given 
by 

P{B Ul <E [a!,bi],...,B Un G [a n ,b n ]) 

dxi--- dx n J^Pui-ui-! (xi-i,Xi) , (Al) 



where uq = 0, xq — 0, and 
p u (a,b) = — 



: exp 



(b-af 
2u 



Eq. ijXTJ can be utilized to compute the expected val- 
ues of moments of standard Brownian motions starting 
at zero. For example, 

E [(B u ) 4 ] = [ -L= exp (~) x 4 dx = Zu 2 , 

where we have used the fact that B u is a Gaussian vari- 
able centered about origin and of variance u, as follows 
from Eq. IjAljl . Therefore, 



E 



/ (B u ) 4 du [ ; [(S u ) 4 ] du = [ 3u 2 du = 1. 
Jo Jo Jo 



(A2) 
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A standard Brownian bridge {£>°, < u < 1} is defined 
as a standard Brownian motion starting at zero that is 
also conditioned to end up at zero at time u = 1: 

{B°,0 < u < 1} = {B u ,0 < u < l|Bi = 0}. 

A standard Brownian bridge can be constructed from a 
standard Brownian motion starting at zero as the differ- 
ence B u — uB\. More precisely, it can be demonstrated 
that 

{Bl,0<u<l} = {B u -uB 1 ,0<u<l}, 

where the symbol = means that the left- and right-hand 
side processes are equal in distribution (have equal finite 
dimensional distributions) and have continuous sampling 
paths with probability one. Moreover, the random vari- 
ables B\ and B® = B u — uB\ are independent. It follows 
that given a Brownian bridge B® and an independent 
standard normal variable z (which plays the role of B\), 
the sum of independent variables B u = B® + uz is equal 
in distribution to a standard Brownian motion starting 
at zero. Thus, 

{B u ,0< u < 1} = {Bl + uz,0 < 1 < u} 

and z = Bi (because B® = 0, by the very definition of 
the Brownian bridge). 

As Simon often emphasizes, 3 Eq. Q presented in the 
introduction is only one of the many equivalent formula- 
tions of the Feynman-Kac formula. Another popular for- 
mulation, which utilizes the full Brownian motion rather 
than the Brownian bridge, will be presented shortly. Let 
ip{x) be an arbitrary square integrable function. From 
Eq. we have 



( 



x e 



-0H\ 



dx'- 



1 



V2? 



exp 



(x' - xf 
2o^~ 



xEexpj-/? J V x + (x' - x)u + aB° u duj tp(x'). 

Performing the substitution x' = x + az, we obtain 

1 



[x e 



-0H\ 



exp (~z 2 /2) 



x exp 



-0 / V 
Jo 



oB\, 



du > tp(x + <rz). 



Notice that the variables z and B®, as they appear in 
the preceding equation, are independent. Moreover, z is 
a Gaussian variable centered in zero and of mean 1. It 
follows that zu + B® is equal in distribution to a Brown- 
ian motion B u starting at zero. In these conditions, the 
Feynman-Kac formula reads 

(x |e"^| i>) = E \ e -Pti v( - x+ ° B ^ du ip(x + aBi)] , (A3) 

where the symbol E denotes the expected value with re- 
spect to the entire Brownian motion B u . 



I conclude this appendix by presenting the statement 
of the Ito-Nisio theorem; 3 ^! theorem that gives an ex- 
plicit construction of a standard Brownian motion over 
the interval [0, 1] random series. 

Theorem 5 (Ito-Nisio) Let {Afc(r)}fc>o be any or- 
thonormal basis in L 2 [0, 1], let 



Afc(u) 



A fc (r)dr, 



and let a :— {ao, a\, . . .} be a sequence of distributed stan- 
dard normal random variables. Then, the random series 
Sfc°=o a k-^k(u) is uniformly convergent almost surely and 
equal in distribution over the interval [0, 1] with a stan- 
dard Brownian motion B u starting at zero. 

To express the Feynman-Kac formula as the expected 
value of a functional of a random series, it is convenient to 
work with those orthonormal basis {\k{T~)}k>a for which 
Aq(t) = 1 only. Then Aq(u) = u and 



A fc (l) 



\k{r)dT 



X k {T)X {T)dT = 



for all k > 1. In these conditions, the Ito-Nisio theorem 
says that 

oo oo 

^ a/cA fc (w) = ^ afcAfc(w) - a Q u = B u - uB x . 



k=0 



The last term in the preceding equation has been dis- 
cussed in a previous paragraph to be equal in distribu- 
tion to a Brownian bridge. It follows that if Ao(r) = 1, 
then 

oo 

S°^^a fe A fe (u), 0<u<l, 

k=l 

equality in distribution that provides an explicit random 
series construction for the standard Brownian bridge. 

In these conditions, if f2 is the set of all sequences d := 
{ai, a,2, ■ ■ ■} and if 



oo 



k=l 



dak 



is the probability measure on f2 associated with 
the sequence of independent random variables a := 
{ai, 02, . . then the Feynman-Kac formula given by 
Eq. reads 



P(x,x';f3) 
p fp (x,x';f3) 

x exp < —/3 



V 



dP[a] 



.(u) + a k A k (u) 



k=l 



du } . (A4) 



Eq. (|A4|) is called the random series representation of the 
Feynman-Kac formula. 10 
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APPENDIX B: THE CONVERGENCE 
CONSTANT FOR THE TRAPEZOIDAL 
TROTTER APPROXIMATION 

The short-time approximation for the trapezoidal 
Trotter path integral method is given by the expression 



-0 



V{x) + V(x') 



This short-time approximation is of the type given by 
Eq. jTSjl , provided that the quadrature technique is speci- 
fied by the two points uq — and u\ = 1, and the weights 
wo = 1/2 and w\ — 1/2, respectively. The approxima- 
tion is independent of the functions {Afe(w);0 < k < q}, 
because the end points of these functions are specified by 
Eq. PJ . We can therefore consider that the functions 
are those for the third order reweighted approximation, 
or one may work with a full random scries representation 
of the Brownian motion of the type 



a u + y^^a k A k (u), 



k=l 

as provided by the Ito-Nisio theorem. It does not make 
any difference. The trapezoidal Trotter approximation is 
just a discrete version of the third order reweighted tech- 
nique discussed in Section IV. A or of the full Feynman- 
Kac formula. 

Using the fact that the trapezoidal quadrature rule 
given above integrates exactly the polynomials 1 and it, 
the reader may argue that all equations specified in Ta- 
blc[I]with n = 1,2,3 are satisfied, except for the following 
(for all, p = 3): 

1) Case je = 1. For the full Brownian motion, one 
computes 



E(M 4 



I [ {B u fdu = [ E(B u ) 4 du = 3 / u 2 du = 1. 
Jo Jo Jo 



The trapezoidal rule produces the different result 



2) Case js = 1 and j\ 
motion, one computes 



1. For the full Brownian 



E(BxM 3 ) = E 



Bx I {B u ydu 
o 



— 31 u du = 1. 
'o 



The trapezoidal rule produces 



3i<o* 



3) Case j'4 = 1 and ji ~ 2. For the full Brownian 
motion, we have 



E[(Bi) 2 Af 2 ] = E 



(B1) 2 / (B u ) 2 du 



2 1 

(2u 2 +u)du =3 + 2' 



The trapezoidal rule produces 
J(2 + l) = 



4) Case j'3 = 2. For the full Brownian motion, we have 
[see Eq. ®] 



Epfi) 2 ] = E 



B u du I = — + — 



1 

12' 



The trapezoidal rule produces 



E 



-1 2 



-(0 + B,) 



With the help of the series given by Eqs. (|2*2"|> and JSJJ, 

we compute 



[p^(x, x';0)- p(x, x'; /})] Mx')dx' = li 3 



4! 1 2 



3! 

2!2! 

(-1) 2 
2! 



1 - 



3, 

1_\ tf_ 

12 J m 



to 
too / 
too 



V {4 \x)^{x) 
V^{x)^ x \x) 
V (2) {x)^ 2 \x) 



V^(x) ^{x)\+0(P 4 ) 



From the equation above, we learn that the trapezoidal 
Trotter path integral technique has convergence order 2. 
Moreover, the convergence operator for the trapezoidal 
Trotter short-time approximation is 



2 48 ImJ [ ' 24 m 



1 fh 2 Yd 



V<»( 3 



._ _ JL v (2) (x) — 
12 \m J dx V dx 



(Bl) 



The above form of Eq. (|B1() emphasizes the Hermiticity 
of the convergence operator. According to Theorem [21 
the following result is expected to hold. 

Theorem 6 The convergence constant for the trape- 
zoidal Trotter path integral method is given by the for- 
mula 

km (n + l) 2 [pl T {x, x'; /?) - p(x, x'; /?)] = 



x')d,8, (B2) 



where the operator T2 is defined by Eq. \hTjj . 

For the purpose of numerical verification, we derive the 
convergence constant for the partition function. Though 
one can work with the full density matrix and employ 
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the Bloch equation whenever necessary, it seems that it 
is more convenient to utilize an eigenfunction expansion 
for the density matrix. Setting x' = x and integrating 
over x in Eq. I|B2|) . we obtain, after several simplifications 
and an integration by parts, 



lim (n+1) 2 [zr{P)-Z(J3)\ 



1 fe2o3 00 



24 m 



fc=0 



2mo 



V (4) (x)ip k (x) 2 - \v {1) {x 



(B3) 



oh 2 

x^ fc (x) 2 + — V^(x) 
mo 



d.j; 



dx. 



Integrating by parts three times, one argues that 

/ V^\x)^{xfdx= — [ V^\x) 
2to J r too J m 



dx 



d 2 

^k{x)—%l} k (x) 



dx, 



whereas, integrating by parts once, we obtain 



o 




100 150 
n 



FIG. 3: The convergence constant for the relative error of the 
partition function of the quartic oscillator, computed for the 
trapezoidal Trotter path integral method. The sequence of 
observed convergence constants c n is seen to converge to the 
theoretical value of c t h ~ 88.35, which is the value predicted 
by Corollary 



too Jr 



Tx Mx) 



dx 



V {1) (x) 



dx 



dx 



ipk{x) 



2tf 
too 



dx. 



Adding the last two equations and simplifying, we get 

f V W(x)Mx) 2 dx+— [ W{x) 
2m J R too J r 



dx 



ipk(x) 



dx = — I V {1) (x) 
too 



d 3 



x i ^k(x)-j^tpk(x) 



-i-ipk{x) 
dx 



<i 

dx 



(B4) 
dx, 



However, by virtue of the Schrbdinger equation, we 
have the equality 



h 2 [ d 3 

d 



dx 



ipk(x) 



dx^ Mx) 



Mx)^{[Ek-V(x)]Mx)} 
dx 



dx 



x {[E k - V{x)]ip k {x)} = -MxfV (1) (x). 
Replacing the last equality in Eq. (|B4(1 . we obtain 



~ / V^(x)Mx) 2 dx+— I V^(x) 
2m Jk to 



d 

dx 



4>k(x) 



dx = 2 



vW(xj\ Mx) 2 dx, 



relation that, upon substitution in Eq. (|B3|) . produces 
the following corollary of Theorem [SJ 



Corollary 1 The convergence constant for the relative 
error of the partition function for the trapezoidal Trotter 
path integral technique is given by the average 



lim [n + 1) 



2 Z ^((3)-Z((3) 



1 h 2 p 3 J R [V^(x)]p(x;f3)dx 



(B5) 



24 to J R p(x;P)dx 
where p(x\ 0) = p(x, x; (3) is the diagonal density matrix. 

Observation. It can be shown that for a multidimen- 
sional system, the convergence constant is given by the 
formula 

ra— >oo Z(p) 

^ 2 P 3 y 1 S R AdiV(x)] 2 P {x-P)dx 
24 £^ m ,i / Kd p(x; 0)dx 

where diV{x) denotes the partial derivative with respect 
to the coordinate i. 

The numerical verification of Corollary ^ is done by 
numerical matrix multiplication for the systems discussed 
in Section V. The theoretical convergence constants 



Cth 



1 h 2 ^J R [V^(x)] 2 p(x-(3)dx 
24 to J R p(x;(3)dx 



can also be computed by numerical matrix multiplication 
(or, more generally, by Monte Carlo integration). The ex- 
perimental values are obtained by numerically studying 
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FIG. 4: Same as in Fig. |3] but for the He cage problem. 



the limit of the sequence 



~ [ + j Z(8) 



As Figs. |3 and 0| show, the agreement between the theo- 
retical and the experimentally observed convergence con- 
stants is excellent for both the quartic oscillator and the 
He cage problem. This agreement is further evidence that 
the statement of Theorem [21 is correct, at least for the 
class of potentials and short-time approximations con- 
sidered in the present article. 
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